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We use ab initio molecular dynamics (AIMD) and the DFT+U method to compute the hydration 
environment of the manganese ion in manganese (II) and manganese (III) porphines (MnP) dispersed 
in liquid water. These are intended as simple models for more complex water soluble porphyrins, 
which have important physiological and electrochemical applications. The manganese ion in Mn(II)P 
exhibits significant out-of-porphine plane displacement, and binds strongly to a single H2O molecule 
in liquid water. The Mn in Mn(III)P is on average coplanar with the porphine plane, and forms a 
stable complex with two H2O molecules. The residence times of these water molecules exceed 15 ps. 
The DFT+U method correctly predicts that water displaces NO from Mn(III)P-NO, but yields an 
ambiguous spin state for the MnP(II)-NO complex. 



I. INTRODUCTION 

Porphyrins dispersed in water, coated on electrodes, 
and self-assembled into nanotubes have important ap- 
plications as sensorsi^^i^ and light-induced water- 
splitting and hydrogen production, 7 Simple manganese 
porphyrins have been used as readily isolated experimen- 
tal models for studying isoelectronic Fe(III) systems." 
Water soluble Mn(II) and Mn(III) porphyrins are par- 
ticularly useful for detection of nitric oxide in cell tis- 
sues, and for distinguishing them from nitroxyl (NO - 
and HNO) compoundsi 2 -^ The relaxation processes of wa- 
ter soluble manganese and iron teira-p-sulfonatcphcnyl 
porphyrins (TSPP) in NMR-compatible time scales 
have also been the subject of significant experimental 
studies,V i u i 12 

Despite their importance and widespread appearance 
in technological and biological settings, theoretical stud- 
ies of transition metal porphyrins in aqueous media have 
been rare. This is in part due to the difficulty of mod- 
cling transition metal systems in general. In liquid wa- 
ter, bare divalent and trivalent first row transition metal 
ions are octahedrally 6-coordinated due to their partially 
filled 3d electron shells. While density functional theory 
(DFT) successfully predicts this feature^ sophisticated 
classical force fields with 3-body terms or quantum me- 
chanics/molecular mechanics methods are needed to re- 
produce such hydration structures, 14 i 15 Transition metal 
porphyrins present an inherently interesting case. The 
molecular framework has an overall -2e charge resulting 
from the nitrogen atoms chelating the metal ions. How 
these chelated ions in porphyrins interact with water will 
be affected by the ion size, spin state, and porphyrin con- 
formations. 

Predicting the spin state of transition metal porphyrins 
presents a considerable theoretical challenge. While DFT 
is formally exact, practical implementations depend on 
the choice of the approximate exchange correlation func- 



tional, not all of which yield the correct spin state in 
transition metal complexes^ This is particularly true 
for first row transition metal porphyrins^ It has been 
shown that the hybrid functionals like B3LYP— yields 
the correct stable high spin ground state for manganese 
(II) porphine (Mn(II)P) and manganese (III) porphine 
(Mnlllljp),! 9 . while non-hybrid functionals like PBE^P. 
lacking a long range Hartree-Fock exchange component, 
predict the incorrect intermediate spin state for Mn(II)P. 
If constrained to the high spin state, PBE still yields 
a poor Mn(II)P molecular geometry , ' 21 However, hy- 
brid functionals are far more expensive to apply than 
non-hybrids when applying periodic boundary conditions 
appropriate to condensed phase systems such as liquid 
water or metal surfaces. This problem has been cir- 
cumventecU 9 - by applying the DFT+U method on the 
Mn 3d electrons.— The screened coulomb ("[/") term 
in DFT+U increases repulsion between electrons in low- 
spin, partially filled d-electron systems. In the ap- 
proach used in Ref. 19, this term was parameterized using 
B3LYP spin splittings, and then the technique was ap- 
plied to model MnP adsorbed on gold electrodes. The 
DFT+U method has also proved accurate for Fe-based 
catalysts.— 

In this paper, we extend this previous work by conduct- 
ing DFT+U based ab initio molecular dynamics (AIMD) 
simulations of Mn(II) and Mn(III) porphines in liquid 
water. Porphines not substituted with ionic side groups 
are mostly water-insoluble; our MnP molecules are in- 
tended as simple models for soluble porphyrins such as 
(MnTSPP)^ ' 10 ! 11 which are too large and computation- 
ally costly to simulate using AIMD £2- The elucidation 
of the hydration environment of the metal ion in por- 
phyrins is a prerequesite for studying the binding of ad- 
ditional ligands to the metal site of water-soluble por- 
phyrins. AIMD hydration studies are potentially useful 
for modeling electrochemical half-reactions involving por- 
phyrins; such calculations have already been performed 



for transition metal ions solvated in water.— Our work 
also suggests that the metal ion hydration structure as- 
sumed in the analysis of NMR relaxation experiments on 
Mn(II)TSPP may need to be re-examined^^ 

Finally, we briefly consider the binding between 
Mn(II)P and Mn(III)P and nitric oxide (NO) in the gas 
phase and in water. Such complexes are classified accord- 
ing to the format {MNO}™^ where n is the number of 
metal d-electrons plus the unpaired electron contributed 
by NO. For n < 6, the M-NO bond should be short 
(~ 1.65 A), and the M-N-0 angle should be ~ 180°^ 
Mn(II) porphyrins are n = 6 complexes, and X-ray struc- 
tures of Mn(II) porphyrins ligated to NO indeed exhibit 
bond lengths and angles appropriate to {MNO} 6 i&2! 
To our knowledge, no X-ray structures of NO-ligated 
Mn(III) porphyrins have been reported. Water soluble 
Mn(III) porphyrins are found not to complex with NO in 
aqueous conditions^ As a result, electrochemical switch- 
ing of the Mn oxidation state is potentially useful for de- 
tecting NO in cell tissues, and in distinguishing NO from 
HNO, which binds to Mn porphyrins in both oxidation 
states. 

Ideally, the same DFT method and/or functional can 
be used to predict all properties (energetics, structures, 
spectral properties) in all aqueous phase applications. 
There have been several reported successes of DFT treat- 
ment of Fe porphyrin-NO complexes ! 28 ' 29 These studies 
have focused on the molecular structures and vibrational 
properties, and not the binding energies between por- 
phyrins and nitric oxide. In this work, we critically exam- 
ine the effect of applying different exchange-correlation 
functionals on the spin state, molecular geometry, and 
binding energies of MnP-NO. We find that the DFT+U 
method successfully predicts that the Mn(III)P-NO com- 
plex dissociates in liquid water environments. How- 
ever, we also conclude that neither DFT+U nor existing, 
widely used DFT functionals enjoy universal succcess in 
modeling nitric oxide ligation. This suggests that sub- 
stantial functional development and refinement are nec- 
essary. 

This paper is organized as follows. Section 2 describes 
the method used. Section 3 discusses the structure of 
water around the manganese ions in porphines. MnP-NO 
binding is briefly described in Sec. 4, and Sec. 5 concludes 
the paper with further discussion. 



II. METHODS 

We conduct ab initio molecular dynamics (AIMD) 
based on the DFT+U2i method implemented^ in the 
VASP code version Our VASP calculations apply 

the PBE exchange correlation functional, the projected- 
augmented wave (PAW) method^ 2 - and the standard 
VASP suite of PAW pseudopotentials. These include the 
Mn pseudopotential with pseudovalent 2p electrons, and 
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FIG. 1: Representative snapshots of (a) & (b) Mn(II)P in 
water; (c) Mn(III)P in water. Except in panel (a), only the 
water molecules close to the Mn ion are depicted. 



the oxygen, carbon, nitrogen, and hydrogen pseudopo- 
tentials with default energy cutoffs of 400, 400, 400, and 
250 eV, respectively. 

The Mn 3d orbitals are augmented with screened col- 
umn ("[/") and exchange ("J") terms of U = 4.2 and 
J = 1.0 eV, respectively. These values have been shown 
to reproduce Mn(II)P high-spin/intermediate spin split- 
tings predicted by the B3LYP functional^ and predict 
the high spin Mn(III) porphyrin spin state observed in 
experiments. For more details, see Ref. 19 and the sup- 
porting information therein. 

A plane wave energy cut-off of 400 eV and a wavefunc- 
tion convergent criterion of 10~ 5 eV at each of the 0.5 ps 
time step in the Born-Oppenheimer AIMD simulations. 
This work focuses on structural properties; hence, deu- 
terium is substituted for all protons so that a larger time 
step can be used. A Nose thermostat with an estimated 
time constant of 20 fs maintains the average tempera- 
ture at T=375 K. With these parameters, the drift in 
total energy is 0.9 Kelvin per picosecond (i.e., ~ 0.004 
mHt/atom/ps) or less. This small energy drift is ab- 
sorbed by the thermostat. Since the underlying PBE ex- 
change correlation functional yields overstructured water 
and slow dynamics at T=300 the elevated temper- 
ature is necessary to obtain reasonable water structure 
and diffusion timescales. (See discussions in the next 
section.) These factors imply that the trajectories do 
not capture the real-time dynamics. Instead, the time 
dependence observed is used to confirm that the system 
has reached equilibrium, and it provides a time scale for 
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thermal fluctuations under these simulation conditions. 
Trajectory lengths of at least 15 ps are used to collect 
statistics after equilibration runs have stabilized the po- 
tential energy. 

The simulation cell is cubic with a linear dimension 
of 13.65 A. The planar porphine molecule is inserted di- 
agonally into the cell. This ensures that there are 5 to 
6 layers of water molecules separating the periodically 
replicated MnP in the direction perpendicular to the por- 
phine ring (Fig. 1). The water density in the cell is de- 
termined using a grand canonical Monte Carlo (GCMC) 
simulation, using the Towhee code^ SPC/E rigid water 
molecule force fields j2£ and porphine force fields appro- 
priate for NiP, 3 §. which, with the exception of the metal 
ion parameters, should be reasonable for porphine atom- 
H2O interactions. The Mn(II)P is held frozen in place 
in the GCMC simulations. It is found that 70 water 
molecules reside in the simulation cell. In the absence of 
Mn(II)P, the 13.65 A 3 cubic cell would contain 85 water 
molecules at 1.0 g/cc density. 

The water content for Mn(II)P-NO is similarly deter- 
mined using GCMC, resulting in 68 H2O molecules in the 
simulation cell. The nitric oxide force field used is sim- 
ilar to that of Ref. 37. NO is essentially a hydrophobic 
molecule that interacts weakly with water. 

All liquid phase pair correlation functions g(r) are com- 
puted with a 0.1 A bin size, except that a 0.05 A bin 
size is used for pure liquid water. 

B3LYP hybrid functional^ calculations of gas phase 
porphine complexes are performed using the program 
Gaussian 03M- We first optimize the geometries using 
the LANL2DZ basis set. Then we switch to the more 
accurate 6-311+G(d,p) basis, and perform geometric op- 
timization for 3 to 5 steps until the magnitude of each 
Cartesian component of the force on each atom is below 
~ 0.05 eV/A. (VASP gas phase calculations utilize a sim- 
ilar criterion for attaining optimal geometry.) We find 
that the PBE Mn(II)P-H20 binding energy computed 
using the 6-311+G(d,p) basis set, and that computed us- 
ing the VASP code with plane wave basis sets, agree to 
within 24 meV. 



III. HYDRATION STRUCTURES OF 
MANGANESE PORPHINES 

First we consider the hydration structure of Mn(II)P 
in water. Simple, electrostatic-plus-Lennard Jones in- 
teractions between water and the porphine molecule are 
applied during the GCMC pre-equilibration simulation. 
These force fields predict that a H2O ligates to the metal 
ion on each side of the porphine plane, but the overall 
metal coordination structure is strongly distorted from 
the ideal octahedral geometry. Applying AIMD on this 
classical force field-generated starting configuration leads 
to substantial structural changes. We find that a 8 ps 




FIG. 2: Pair correlation functions g(r) between the Mn ion 
and the oxygen sites of H 2 for (a) Mn(II)P; (b) Mn(III)P. 
Inset: goo(r) for liquid water. 



AIMD trajectory is needed to equilibrate the potential 
energy of the system. AIMD simulation of Mn(III)P in 
water is started from an equilibrated Mn(II)P configu- 
ration, and the potential energy of the system also con- 
verges to a plateau value within 8 ps. These long equili- 
bration times help ensure that the statistics collected in 
the subsequent 15 ps production run are independent of 
initial conditions. 

Figure la depicts a equilibrated hydration structure 
around Mn(II)P and confirms that there are several lay- 
ers of water molecules buffering the Mn ion from its pe- 
riodic image. This suggests that the hydration environ- 
ment of the Mn ion adequately represents bulk water 
boundary conditions. 

The representative configuration depicted in Fig. lb 
focuses on H2O close to the Mn ion. Unlike the classical 
force field results, the Mn ion in Mn(II)P strongly binds 
to one rather than two H2O molecules. As discussed be- 
low, this can be related to the significant out of porphine 
plane motion exhibited by the Mn(II) ion. Figure lc de- 
picts Mn(III)P in water. Mn(III)P is on average coplanar 
with the porphine plane, and the Mn binds to two wa- 
ter molecules. These snapshots encapsulate this section's 
main results, which will be analyzed in some detail below. 

The Mn(II)-H 2 pair correlation function <7Mn-o w M 
is plotted in Fig. 2a, and it confirms the structure shown 
in Fig. lb. The first peak in the Mn-O w g(r) integrates 
to 1.0 from r = to r = 3.0 A, which is the first g(r) 
minimum. The sharp peak and deep minimum clearly 
show that this water molecule is strongly bound to the 
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Mn(II) ion. They are reminiscent of divalent cation g(r) 
in liquid water i 15 i 39 In fact, over the entire 15 ps pro- 
duction run trajectory, the same water molecule remains 
bound to the Mn ion, although during the 8 ps equili- 
bration trajectory, exchange of Mn first hydration shell 
water molecules with the bulk water region has occurred. 
Experimentally, the time scale for proton exchange that 
involves H2O bound to Mn(II) is known to be on the 
order 10" 8 to 10" 7 si2 

The inset to Fig. 2 depicts the pure water goo{ r ), com- 
puted using the PBE functional, 32 water molecules at 
1.00 g/cc light water density, and with the thermostat 
set at T=375 K. It confirms that PBE water at T=375 K 
exhibits less overstructuring compared to PBE water at 
T=300 K. For example, the goo{r) first peak exhibits a 
peak height of 2.9 density units at T=375 K, which is 
in reasonable agreement with the value of 3.1 reported 
by Sit et al.£^ The small discrepancy can be attributed 
to statistical uncertainties and/or the use of different bin 
sizes. In contrast, at T=300 K, this peak height exceeds 
3.4 density units^ which is considerable larger than the 
experimental value. 41 Our choice of T=375 K is thus 
reasonable for MnP hydration studies. While we have 
imposed a water density estimated using classical force 
fields and the GCMC method, we note that ion hydra- 
tion structure appears insensitive to significant variations 
(~ 7%) in water density^ 

The stability of the single- H2O coordination structure 
of Mn(II)P is closely related to the significant Mn out- 
of-plane motion. Figure 3a depicts the evolution of this 
displacement along the AIMD trajectory. Here, at each 
time step, we perform a least square-fit to generate a 
plane, ax + by + cx = 1, through the 4 nitrogen atoms, 
and <5mii is obtained as the smallest distance between the 
Mn atom and this plane. On average, <5mii is 0.48 A for 
Mn(II)P dispersed in water. 

This large displacement can be related to the extremely 
flat out-of-plane Mn potential energy surface predicted 
by the DFT+U method for high spin Mn(II)P in the gas 
phase (Fig. 3b). Thus, upon binding to a ligand like 
a H2O, the Mn ion is readily displaced out of the por- 
phine plane. We find that this displacement is already 
0.225 A for the Mn(II)P-H20 complex in the gas phase, 
and this geometric feature explains the inability of the 
Mn ion to bind to a second H2O in a liquid water enviro- 
ment. Note that the out-of-plane displacement of the Mn 
ion is accompanied by only a small amount out-of-plane 
deformation of the porphine ring (approximately 0.1 A 
total deformation as determined by Normal Coordinate 
Structure Decomposition; see Table SI). 

The out-of-plane Mn motion exhibits a flat potential 
energy surface because the high spin Mn(II) ion has a 
large effective size. The Mn-N distance in high spin 
Mn(II)P is predicted to be 2.09 A, in good agreement 
with experimentsi 45 i Thus, energetically, it is relatively 
unfavorable for Mn to fit in the center of the ion-chelating 



site of the porphine molecule. In contrast, the Mn(III) 
ion is smaller, with an Mn-N distance of 2.045 A in 
MnClPfi^ and in isolated Mn(III)P strongly prefers to 
be coplanar with the porphine plane (Fig. 3c). In liq- 
uid water, this preference persists on average (Fig. 3a). 
This geometry allows effective binding to two ligands. 
Indeed the sharply peaked Mn- water first peak in g(r) 
integrates to 2.0, indicating two H2O are strongly lig- 
ated to the Mn ion. We do not observe exchange of the 
two H2O molecules bound to the Mn(III) ion with outer 
shell water during the 15 ps trajectory. 

A search of the Cambridge Crystallographic Data Cen- 
ter (CCDC) yields 15 crystal structures of Mn porphyrins 
with one or more coordinated water molecules. The man- 
ganese ion exhibits the +3 oxidation state in all of them. 
The experimentally determined Mn-0 distances in the 
Mn(III) porphyrins (~ 2.2 A) correspond closely to the 
first peaks in the pair correlation functions in Fig. 2. 
The average calculated Mn-N (porphine) distance in the 
high-spin (m = 2) Mn(III)P(H2 0)2 energy minimized 
using the DFT+U is 2.028 A. This again agrees well 
with the Mn-N(porphyrin) distances observed in Mn(III) 
porphyrin complexes with two axial water ligands (e.g., 
Mn(III)TPP(H 2 0)2(BPh)4, where the Mn-N distance is 
2.011 A) 41 The porphine ring again shows very little 
non-planar distortion in either the DFT+U gas phase 
structure of Mn(III)P(H 2 0)2 (Table S2) or the crystal 
structure of Mn(III)TPP(H 2 0)2(Bph)4 (Table S3). 

Gas phase B3LYP and DFT+U calculations of the 
MnP-(H 2 0) 2 complex (Fig. 4, Table QJ provide more 
compelling evidence that Mn in Mn(II)P coordinates to 
only one water molecule. Both these methods predict 
that there are two nearly isoenergetic configurations for 
the high spin Mn(II)P-(H 2 0) 2 cluster (Fig. 4). The 
bidentate Mn(II) configuration, with the Mn ion resid- 
ing in the porphine plane (Fig. 4a), is only slightly more 
stable than the one where Mn(II) is ~ 0.4 A out-of- 
plane and one H2O is bound to Mn and the other co- 
ordinated to two nitrogen atoms (Fig. 4b). DFT+U 
predicts that the bidentate configuration is more stable 
by only 15 meV, while B3LYP predicts that they are 
equally stable. In contrast, with the system in either the 
intermediate or the high spin state, PBE predicts that 
the structure depicted in Fig. 4b is unstable and sponta- 
neously relaxes to that shown in Fig. 4a. This emphasizes 
that there are fundamental differences between DFT+U 
and PBE predictions for gas phase properties. These dif- 
ferences carry over into liquid phase AIMD simulations 
(Fig. la) favoring configurations similar to Fig. 4b4 8 - 

Because of the slow exchange rate of Mn first hydra- 
tion shell H2O, it is important to check whether the 
Mn(II)P hydration structure we observe depends on the 
initial configurations. Thus, we start with a equilibrated 
Mn(III)P configuration, with two water molecules lig- 
ated to the Mn ion, and initiate a Mn(II)P AIMD sim- 
ulation by adding one more electron to the simulation 
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FIG. 3: (a) Out-of-porphine plane displacement of Mn ion, 
<5 M n, over the AIMD trajectories for Mn(II)P and Mn(III)P. 
(b) & (c) Energies as functions of Mn out-of-plane displace- 
ment in the gas phase for Mn(II)P and Mn(III)P, respectively 
The porphine ring is held fixed as the Mn is displaced. 



cell. Within 1 ps, one of the first hydration H2O de- 
taches itself, the Mn ion exhibits strong out-of-plane dis- 
placement, and the single H2O hydrated, 5-coordinated 
Mn(II) configuration is recovered. Thus, this Mn(II)P 
hydration structure is robust. 

The predicted hydration structure of Mn(II)P might 
possibly resolve anomalies associated with paramagnetic 
relaxation behavior of Mn(II)TPPS in water . 10 ' 11 One 
model proposed in the NMR literature assumes a 6- 
coordinated geometry, where 2 axial H2O molecules bind 
to the manganese ion in Mn(II)TPPS. 11 In contrast, 
AIMD predicts only one H2O in the hydration shell and a 
breaking of the inversion symmetry about the porphine 
plane, which now exhibits a C^ v symmetry. This sug- 
gests that the underlying assumptions of Ref . 1 1 regard- 
ing coordination to water may need to be re-examined. 
While paramagnetic relaxation behavior is complex, we 
propose that using a model with the Mn ion ligated to 
the AIMD predicted number of water molecules might 
improve agreement between experiments and the spin- 
Hamiltonian results used to inteprete NMR data ] 10 ' 11 

Figure 5 examines the hydration environment of the 
four nitrogen atoms surrounding the Mn ion. Integrat- 
ing the <7N-o w ( r ) to its first minimum (3.8 and 3.9 A 
for Mn(II)P and Mn(III)P respectively), we find that 
there are 1.76 and 2.83 H2O in the hydration shell of 
each nitrogen atom. These spatial correlations partly 
originate from the one (two) H2O strongly coordinated 
to the Mn(II) (Mn(III)) ion in the vicinity of the nitro- 
gen atoms, and not from water-nitrogen hydrogen bond- 
ing. A N-H hydrogen bond typically exhibits a <7n-h (t) 
peak at r ~ 1.8 A, and a first minimum at r ~ 2.5 A. 49 



(a) 




(b) 



FIG. 4: Two configurations of Mn(II)P-(H 2 0) 2 . DFT+U and 
B3LYp both predict that they are similar in energy, while 
PBE predicts that the structure in (b) is unstable. 



Here, gN-n^(r) only exhibits small shoulders. Integrat- 
ing <7n-h„M to r = 2.5 A yields only 0.10 and 0.02 
H2O in the first coordination shell per nitrogen atom for 
Mn(II)P and Mn(III)P, respectively. We conclude that 
water interacts weakly with the nitrogen atoms. 

Note that our AIMD trajectories impose a high spin 
configuration, known to be favorable in the gas phase. 19 
We have tested snapshots of the Mn(II)P trajectory to 
see if Mn(II)P in water favors spin states at variance from 
those predicted for the gas phase ground state. We invari- 
ably observe that the high spin configuration is favorable 
in these DFT+U generated AIMD snapshots when using 
the DFT+U technique, in agreement with experiments 
on MnTSPP in water . 10 ' 11 The intermediate spin state 
is favored in these snapshots when the PBE functional 
is used. Thus, dispersing Mn(II)P into liquid water does 
not alter the most stable spin state; for each functional 
used, the gas phase and aqueous phase spin state predic- 
tions are in agreement. This is expected because water 
is not a strong ligand. 



IV. INTERACTIONS WITH NITRIC OXIDE 

In this section, we examine the interaction of Mn(II)P 
and Mn(III)P with NO. Ideally, the DFT+U method, 
which predicts the correct s = 5/2 spin state for Mn(II)P, 
can be used to treat MnP-NO binding in aqueous en- 
vironments. We will show that DFT+U appears to 
yield reasonable binding energies. Neither DFT+U nor 
B3LYP necessarily predicts the correct spin states for 
the MnP-NO complexes. However, since the energies of 
different spin states are so similar, our conclusions con- 
cerning binding energies appear to be robust. 

We first consider gas phase predictions. Table[T]depicts 
the PBE, B3LYP, and DFT+U results based on the PBE 
functional and U = 4.2 eV, J = 1.0 cV. Here we use the 
notation m to denote the magnetic moment of the entire 
complex in the simulation cell. 
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TABLE I: Binding energies between MnP and NO or H2O, 
in eV.— PBE and DFT+U results utilizes planar wave ba- 
sis while B3LYP calculations are performed using the 6- 
311+(d,p) basis, m is the total magnetic polarization in the 
complex. In some cases only the results for the most sta- 
ble spin state predicted for each DFT functional are listed. 
'Negligible binding energy predicted. 



A. Mn(II)P-NO 

PBE correctly predicts a low spin ground state for 
Mn(II)P-NO, with a 1.617 A Mn-N (nitric oxide) bond 
length, and a Mn-N-0 angle of 173.6°. The Mn ion 
is displaced 0.342 A out of the porphine plane, and 
the porphine ring shows little evidence ofo non-planar 
deformation (~ O.lA, see Table S4). These predic- 
tions are in good agreement with the X-ray struc- 
ture for nitrosyl(5,10,15,20-tetratolylporphinato) man- 
ganese (Mn(TTP)NO), which reveals a 1.641 A NI- 
NO bond length, a 177.8° M-N-0 angle, and a 
0.337 A Mn out-of-plane displacement The four 
other Mn(II)tetraarylporphyrins in the CCDC are 6- 
coordinated but show similar Mn-NO distances (1.645- 
1.680 A). The average Mn-N(porphyrin) distance in the 
PBE structure of Mn(II)P-NO is 2.022 A, which com- 
pares favorably with the crystal structure distance of 
2.004 A. 26 The crystal structure of Mn(II)TPP-NO also 
shows a modest amount of non-planar deformation (~ 
1.2 A, Table S5) which may in part explain the slightly 
shorter Mn-N (porphyrin distance. PBE predicts a large 
(2.03 eV) MnP-NO binding energy (Table U). 

DFT+U and B3LYP results are qualitatively and even 
quantiatively similar to each other but disagree with the 
experimental data. They predict, incorrectly, that the 
m = 1, intermediate spin configuration is more stable 
than the low spin state, by 0.20 and 0.32 eV, respec- 



r(A) 

FIG. 5: Pair correlation functions, <7n-o w (?") (solid lines) and 
ffN-H.(f) (dashed lines), between the nitrogen atoms and the 
oxygen/hydrogen sites of H2O for (a) Mn(II)P; (b) Mn(III)P. 



tively. For this ground state, the Mn-Nwo distance is 
1.91 (1.89) A the Mn-N-0 angle is 142.2 ( 144.2)°, the 
out-of-plane Mn displacement is 0.277 (0.323) A, and the 
MnP-NO binding energy is 0.93 (0.73) eV. DFT+U fur- 
ther predicts that the low-spin configuration still exhibits 
a considerable binding energy of 0.61 eV, with a 1.766 A 
Mn-NO bond length and 174° M-N-0 bond angle, re- 
spectively. 

We note that Gaussian B3LYP calculations reveal 
some spin-contamination for the m = 1 intermedi- 
ate spin state. The VASP package does not calculate 
spin contaminations, but it is likely that VASP-based 
DFT+U calculations also contain a similar degree of spin- 
contamination for this spin state. We emphasize that 
only the NO-ligated species exhibit such a contamina- 
tion. Since both the energies of the m = and m = 1 
spin states agree to within 0.3 eV, the order of magnitude 
of our binding energies will likely not be strongly affected 
by the ambiguity regarding the spin state. 



B. Mn(III)P-NO 

To our knowledge, no X-ray structures for NO-ligated 
Mn(III) porphyrins exist. Such hypothetical systems 
would be {M-NO} 5 complexes, which in principle should 
exhibit a low-spin (m = 1/2), short Mii-Nno, linear Mn- 
N-0 geometry. 26 PBE predicts such a low spin structure, 
but also yields an anomalously large MnP-NO binding 
energy of 1.30 eV (Table 12 see below). This structure 
has an out-of-plane Mn displacement of 0.36 A, which is 
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(a) (b) 



species 


expt. 


PBE B3LYP DFT+U 


Mn(II)P 


s=5/2 


no 


yes yes 


Mn(III)P 


s=2 


yes 


yes yes 


Mn(II)P-NO 


s=0 


yes 


no no 


Mn(III)P-NO 


instability in water 


no 


NA yes 



TABLE II: Brief summary of the successes and failures of 
different DFT methods. Gas phase DFT and aqueous phase 
AIMD simulations that apply the same functional predict the 
same stable spin states in all cases considered. 



FIG. 6: Representative snapshots of (a) initial and (b) final 
configurations of the Mn(III)P-NO complex in water after 
0.1 ps. Only H2O molecules close to the Mn ion are depicted. 



comparable to that seen for low spin Mn(II)P-NO. 

Both DFT+U and B3LYP predict negligible binding 
energies for NO with low spin Mn(III)P. The m = 3/2 
intermediate spin configuration, on the other hand, yields 
a small binding energy of 0.40 (0.32) eV and a small out- 
of-plane Mn displacement of 0.16 (0.10) A. While no ex- 
perimental data for binding energy is available, DFT+U 
predicts that the Mn(III)P-H 2 complex is more favor- 
able than the Mn(III)P-NO complex, while the reverse is 
true for PBE. This suggests that H2O will displace NO 
in a DFT+U based AIMD simulation, consistent with 
experiments on water-soluble Mn(III) porphyrins'^ In 
contrast, because of the strong binding energy with NO, 
AIMD simulations based on the PBE exchange correla- 
tion functional will likely not lead to Mn(III)P-NO dis- 
sociation in water 1^ 

Gaussian calculations again reveal some spin- 
contamination for both the m = 1/2 and m = 3/2 
states. In contrast, no spin-contamination is observed 
for either Mn(II)P or Mn(III)P-ClJ£ This underscores 
that MnP-NO systems are difficult to model using DFT 
methods. 



C. MnP-NO in water 

Next, we investigate the stability of MnP-NO com- 
plexes in liquid water by conducting AIMD simulations. 
We pre-equilibrate hydrated Mn(II)P-NO configurations 
from a GCMC run as before, holding MnP-NO rigid in 
the DFT+U optimized, intermediate spin configuration, 
and then conduct AIMD simulations for 3 ps. Unlike for 
Mn(II)P ligated to H2O, the potential energy reaches a 
plateau on a sub-picosecond timescale. We find that the 
Mn ion in Mn(II)P-NO complex is stable in water, and 
does not ligate to any additional H 2 molecules. 

Next, we consider Mn(III)P-NO in water. Starting 
from an equilibrated aqueous Mn(II)P-NO configuration 
(Fig. Id) , we remove an electron from the simulation cell, 
and restart the DFT+U based AIMD trajectory with the 



spin polarization fixed at m — 3/2. Within 0.1 ps, the 
NO diffuses away from the Mn(III) ion, replaced by an 
H2O from the other side of the porphine plane (Fig. le; 
Fig. 6). Thus, DFT+U predicts that Mn(III)P fails to 
bind to NO in water, in agreement with experiments. 2 

To examine PBE predictions for Mn(III)P-NO, we use 
a DFT+U predicted Mn(III)P configuration in water, 
manually reposition the NO molecule to yield a 180° Mn- 
N-0 angle and a 1.8 A Mn-N bond length, impose a low 
spin (to = 1/2) spin state, and restart the AIMD simula- 
tion using the PBE functional. We find that Mn(III)P- 
NO complex is stable in water over the course of a 5 ps 
trajectory. While these AIMD runs are short, their qual- 
itative results are completely consistent with the rela- 
tive stabiliy of Mn(III)P-NO and Mn(III)P-H 2 based on 
gas phase binding energy considerations (Table HI . They 
show that PBE predictions are in apparent disagreement 
with experiments.— 



D. Summary of MnP-NO interactions 

The apparent successes and failures of PBE and 
DFT+U methods in the gas and liquid phases are sum- 
marized in Table [ill Recall that, in the absence of 
the NO ligand, PBE predicts an incorrect, intermedi- 
ate Mn(II)P (s=3/2) spin state while both B3LYP and 
DFT+U parameterized with B3LYP spin splittings pre- 
dict the correct high spin (s=5/2) state, 19 All three meth- 
ods predict the correct spin state for Mn(III)P (high spin, 
s = 2). The DFT+U method correctly predicts that the 
Mn(III)P-NO is unstable in liquid water. However, the 
spin state of the predicted stable species in the gas phase 
is inconsistent with the typical low-spin {M-NO} 5 struc- 
tural parameters. PBE predicts that Mn(III)P-NO is a 
strongly bound complex, in apparent disagreement with 
aqueous phase experiments. For the Mn(II)P-NO struc- 
ture, PBE predicts the experimental structure and spin 
state, while B3LYP and DFT+U do not. 

Both the hybrid functional B3LYP and the DFT+U 
method favor the high-spin states of first row transi- 
tion metal ion, by increasing the exchange interaction or 
screened couloumb interaction among 3c? electrons. To 
a first approximation, when Mn and NO are far apart, 
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enforcing m — on the Mn(II)P-NO complex requires 
equal and opposite local s = 1/2 spin moments on Mn 
and on the N atom of the NO molecule. We conjecture 
that DFT+U cannot reproduce the stability of m = 
Mn(II)P-NO because Mn(II)P strongly favors the high 
spin (m = 5/2) state, while the low spin (m = 1/2) state 
is highly unfavorable. In contrast, in PBE calculations 
the m = 1/2 state is not extremely unfavorable in energy 
compared to m = 3/2, the stable Mn(II)P spin state er- 
roneously predicted by PBE. Therefore, at least when 
Mn(II)P and NO are far apart, PBE predictions more 
readily accommodates a m = spin state. 

In summary, state-of-the-art DFT methods are as yet 
unable to capture the delicate balance of competing ef- 
fects and correlations that determine the stable spin state 
of Mn(II)P, Mn(III)P, and their nitric oxide complexes. 
NO ligation seems to be especially challenging due to ir- 
electron back bonding^^ The DFT+U method, with 
its most important parameter U fitted to the Mn(II)P 
high spin-intermediate spin splitting, should be further 
improved to take into account such 7r-back bonding. 53 
These improvements may allow the more accurate pre- 
diction of binding energies, which have not been ad- 
dressed in recent DFT work on porphyrin-nitric oxide 
complexes j 28 ' 29 

V. CONCLUSIONS 

In this work, we have conducted ab initio molecular 
dynamics simulations of the hydration environment of 
the Mn ion in Mn(II)P and Mn(III)P dispersed in liquid 
water. These are intended as simple models of water- 
soluble manganese porphyrins. The DFT+U technique, 
parameterized with B3LYP spin splittings, successfully 
predicts that high spin Mn(II)P and Mn(III)P are stable 
in liquid water. Thus, this technique enables efficient 
molecular dynamics simulation of MnP with the correct 
spin state in a condensed phase environment. In contrast, 
the B3LYP functional, with its long range exchange, is 
costly to apply when using periodic boundary conditions. 

The Mn ion in Mn(II)P is predicted to displaced off- 
center out of the porphine plane by an average of 0.48 A, 
and is ligated to a single H2O molecule. The Mn ion 
in Mn(III)P is on average co-planar within the porphine 
plane and binds strongly to 2 H2O molecules. Water 
molecules ligated to the Mn ion exhibit residence times 
of more than 15 ps. Water only interacts weakly with 
the nitrogen atoms on the porphine ring chelating the 
Mn ion. These predicted hydration structures might be 
potentially useful for improved analysis of NMR relax- 



ation datai 10 ! 11 

The application of DFT techniques to examine MnP- 
NO binding in both gas and aqueous phases proves less 
successful than for MnP dispersed in water. Ideally, 
the DFT+U technique, fit to B3LYP spin splittings for 
Mn(II)P and found to be successful for unligated MnP, 
should be applicable to all ligated complexes as well. In- 
deed, we find that DFT+U apparently succeeds in de- 
scribing the binding energy of the Mn(III)P-NO complex. 
It predicts that the Mn(III)P-NO binding energy is sub- 
stantially smaller than that for Mn(III)P-H 2 0. AIMD 
based on the DFT+U technique predicts a rapid dis- 
sociation of the Mn(III)P-NO in water, consistent with 
the known instability of Mn(III)P-NO in aqueous phase 
experiments^ 3 - Both DFT+U and B3LYP predict sig- 
nificantly larger Mn(II)P-NO binding energies than for 
Mn(ni)P-NO. This differential ability of Mn(II) and 
Mn(III) porphyrins to ligate to NO is crucial for detect- 
ing NO in water and distinguishing NO from HNO or 
NO~ via electrochemical means. In contrast, PBE pre- 
dicts large (> 0.8 eV) Mn-NO binding energies for both 
Mn(II)P-NO and Mn(III)P-NO, which would suggest, er- 
roneously, that switching the oxidation state of MnP can- 
not trap and then release NO. 

We emphasize that B3LYP and DFT+U predictions 
qualiltatively agree with each other regarding Mn(II)P- 
NO binding, despite the fact that the DFT+U technique 
used herein is parameterized using B3LYP spin splitting 
for the bare Mn(II)P molecule. However, both these tech- 
niques predict incorrect spin states for the Mn(II)P-NO 
complex. On the other hand, the PBE functional pre- 
dicts a Mn(II)P-NO low spin spin state and molecular 
geometry that are in good agreement with experiments. 

While it is clear that more fundamental improvements 
to DFT methods need to be made, AIMD simulations 
based on the DFT+U technique successfully predicts 
many experimentally observed ligation features of Mn 
porphyrins, including the observation that complexes 
formed between Mn(III) porphyrins and nitric oxide are 
unstable in liquid water. 
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